function [x, U] = SquareRootMethod(A, f)
[R, p] = chol(A);
N = size(A, 1);
S = zeros(N);
D = eye(N);
if (p == 0)
    %Positively defined matrix
    for i = 1 : N
        
        sum = 0;
        for p = 1 : i - 1
            sum = sum + S(p, i)^2;
        end;
        S(i, i) = sqrt(A(i, i) - sum);
        
        for j = i + 1 : N
            sum = 0;
            for p = 1 : i - 1
                sum = sum + S(p, i) * S(p, j);
            end;
            S(i, j) = (A(i, j) - sum) / S(i, i);
        end;
        
    end;
else
    %Negatively defined matrix
    for i = 1 : N
        
        sum = 0;
        for p = 1 : i - 1
            sum = sum + S(p, i)^2 * D(p, p);
        end;
        D(i, i) = sign(A(i, i) - sum);
        S(i, i) = sqrt(abs(A(i, i) - sum)); 
        
        for j = i + 1 : N
            sum = 0;
            for p = 1 : i - 1
                sum = sum + S(p, i) * D(p, p) * S(p, j);
            end;
            S(i, j) = (A(i, j) - sum) / (D(i, i) * S(i, i));
        end;       
    end;
end;
%Got matrix S
B = S' * D;

y = zeros(N, 1);

for  k = 1 : N
    sum = 0;
    for s = 1 : k - 1
       sum = sum + B(k, s) * y(s); 
    end
    y(k) = (f(k) - sum) / B(k, k);
end;

x = zeros(N, 1);

k = N;
while (k >= 1)
    sum = 0;
    for p = k + 1 : N
        sum = sum + S(k, p) * x(p);
    end;
    x(k) = (y(k) - sum) / S(k, k);
    k = k - 1;
end;


dlmwrite('lab4_report.txt', S, '-append', 'precision', '%.4f', 'delimiter', '\t');
dlmwrite('lab4_report.txt', D, '-append', 'precision', '%.4f', 'delimiter', '\t', 'roffset', 1);
dlmwrite('lab4_report.txt', y, '-append', 'precision', '%.4f', 'delimiter', '\t', 'roffset', 1);

U = A * x - f;

end

